Vortex quantum creation and winding number scaling in a quenched spinor Bose gas 
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Motivated by a recent experiment, we study non-equilibrium quantum phenomena taking place 
in the quench of a spinor Bose-Einstein condensate through the zero-temperature phase transition 
separating the polar paramagnetic and planar ferromagnetic phases. We derive the typical spin 
domain structure (correlations of the effective magnetization) created by the quench arising due to 
spin-mode quantum fluctuations, and establish a sample-size scaling law for the creation of spin 
vortices, which are topological defects in the transverse magnetization. 
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A rapid sweep through a symmetry-breaking zero- 
temperature phase transition entails fascinating non- 
equilibrium quantum phenomena. The quench gener- 
ates causally disconnected spatial domains of different 
order parameter values (corresponding to positions in the 
new ground-state manifold), which may ultimately form 
topological defects in a quantum version of the Kibble- 
Zurek mechanism [3, 0, 0] . Apart from its relevance in 
condensed-matter theory (e.g., in superconductors, these 
defects affect measurable transport properties like con- 
ductivity and susceptibility) and potentially cosmology 
[2], this phenomenon is very interesting from a funda- 
mental point of view, since the defect distribution di- 
rectly originates from the initial quantum fluctuations 
and thus maps out their properties (such as the ground- 
state entanglement). 

Because of their comparably long response times, di- 
lute atomic gases provide a unique opportunity for ex- 
ploring these fundamental quantum many-body phenom- 
ena far from equilibrium in the laboratory. In the follow- 
ing, we derive the spectrum of fluctuations induced by 
a rapid quench to the (planar) ferromagnetic state of a 
spinor Bose gas, and establish a general sample-size de- 
pendent scaling law for the resulting variance of the net 
number of spin vortices (i.e., of the winding number). Be- 
sides the total defect density, the winding number within 
a given area determines the spectrum of the created mag- 
netization fluctuations. Thus dilute Bose gases serve as 
a useful and experimentally accessible toy model for the 
general case and may allow us to extract universal prop- 
erties of such dynamical phase transitions. 

Spinor Bose-Einstein condensates are created by trap- 
ping different hyperfine states of a particular atomic 
species by optical means [3], which enables, inter alia, 
the investigation of coherent spin-exchange dynamics [f|, 
and the formation of spin domains by a dynamical insta- 
bility Non-equilibrium phenomena in a spinor Bose 
gas were realized in a recent experiment , where an ini- 
tially paramagnetic state was rapidly quenched through 
a quantum phase transition to a final ferromagnetic state 
H . Such a quench results in spin vortices, (imaged in situ 



by phase contrast techniques [9j), which are topological 
defects in the magnetization with a paramagnetic core 
[ToL Hi]]. The dilute spinor Bose gas can be described in 
terms of multi-component field operators ip a , whose dy- 
namics is governed by the Hamiltonian density (ft = 1) 

+ Y Fab ' Fcd ftatilMd ~ idtd* ' (!) 

where m denotes the mass of the atoms. The vector 
F a b contains the spin matrices [13] and determines the 
effective magnetization F = F ab t[)li/; b / g, where a sum- 
mation over component indices a,b,c,d is implied and 
q = (x^ltpa) is the total (conserved) density. For spin-one 
systems, the sum runs over a — 0,±1 or, alternatively, 
over a — x,y,z with ipo = tp z and tp± — (ip x ± j^)/v2- 
The coupling constants co = 47r(ao + 2a2)/(3m) and 
C2 = 47r(<22 — ao)/(3m) are determined by scattering 
lengths as for the scattering channel with total spin S. 
For Co ^> I C2 1 , density fluctuations are energetically sup- 
pressed in comparison with the spin modes. Since we are 
interested in the effects of the phase transition (which is, 
strictly speaking, only well-defined for infinite systems) 
and not in the impact of inhomogeneities, we assume 
g const., and hence omit the scalar trapping potential 
Vt ra p. Finally, q denotes the strength of the quadratic 
(one-particle) Zeeman shift qF 2 [7j, where an external 
magnetic field is oriented along the z direction. (The lin- 
ear Zeeman shift can be eliminated by transforming to 
the co-rotating frame, because the Larmor precession is 
much faster than all other frequency scales [id].) 

Assuming C2 < 0, the system is ferromagnetic for van- 
ishing Zeeman shift q = 0, i.e., the magnetization F as- 
sumes a maximum value (F) 2 = 1 in some given but arbi- 
trary direction (broken-symmetry phase). For small but 
non-zero q, the magnetization (F) lies in the x, y-plane 
due to the term qF 2 but is still maximal, (F) 2 = 1, For 
large q, however, the ground state corresponds to a con- 
densate in the ipo — ip z component only and symmetry 
is restored, i.e., the magnetization (F) vanishes ("para- 
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FIG. 1: Inversion behavior of the effective potential from (f2|). 
After crossing the critical point, the direction n of the effective 
magnetization F in the plane, at any given point in space, gets 
frozen at a particular angle while its modulus F continues 
to grow exponentially in time ("rolls down" the potential hill). 



magnetic" phase). Hence, by rapidly lowering q, we may 
quench the system from an initially paramagnetic to an 
(effectively two-dimensional) ferromagnetic phase. 

In order to study the evolution of the quantum field 
during the quench, we employ a number-conserving 

mean-field ansatz 4>z = (V'co + S^p z )A/\/N , which 
is adapted to the initial paramagnetic phase but can 
be extrapolated for some finite time after the quench 
as long as the quantum fluctuations are small 
enough 8ip x ,Stfj y ,8ip z <§; tp co . This ansatz allows for a 
linearized expression for the magnetization (which as- 
sumes in Cartesian coordinates the simple form F a — 
—i^abci > \i ) c/Q) and to derive an effective mean-field 
Hamiltonian for transverse magnetization F from the full 
Hamiltonian Eq. ([!]), 
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with fl = (ipcotpl + iPco^y , -V'co^J - ip* 4> x )/2 being the 
canonical momentum operator for spin-one bosons. The 
experiment [?J is described by a two-dimensional trans- 
verse magnetization F — (F xl F y ) obeying the 50(2)- 
symmetry of rotations around the z-axis in effectively 
two spatial dimensions [V = (d x ,d y )], permitting topo- 
logical defects in the form of spin vortices. However, 
the above expressions for the Hamiltonian are completely 
analogous for a 50(iV)-symmetry of the order parameter 
(F) in N > 2 spatial dimensions, where the analogous 
topological defects are generically called "hedgehogs", see, 
e.g., EES HQ- Therefore, we shall discuss the gen- 
eral SO(N) situation in the following and return to the 
experimentally realized example N = 2 of 0] at the end. 

The simple expression ([2j) allows us to directly read 
off the critical value q CI = 2 1 C2 1 £> of the phase transition. 
After a quench from q m > q cr to q ou t < q a , and us- 
ing a normal-mode expansion, we obtain exponentially 
growing modes corresponding to imaginary frequencies 
lo\ = V ( e k + q)( e k + q + 2c 2 £o), where e k = fc 2 /(2m), 
cf. Fig.[TJ There are two regimes: For q CT > q out > q C i/2, 
the exponential growth rate ^(wjQ of the normal modes 
increases with their wavelength, but for q ou t < g C r/2, the 



growth rate assumes its maximum at a given wavenumber 
/; max = \J m(q cr — 2g out ). In view of the experimental pa- 
rameters [7||, we shall focus on the second case. Note that 
due to the emergence of this preferred length scale, which 
is independent of the quench time, the present situation 
is somewhat different from the "standard" Kibble-Zurek 
scenario [l|, 0] , which is based on the interplay between 
quench time, relaxation rate, and correlation length. 

Assuming the Bose gas to be sufficiently dilute, we 
may extrapolate our linearized (mean-field) description 
([2]) to intermediate times t which contain many e-foldings 
^(u^ )t ^ 1 of the exponentially growing modes [1J|, 
whereas the magnetization is still small enough, f 2 <l, 
to stay in the linear regime. After a normal-mode expan- 
sion, the expectation value (F a (r)Fb(r')) contains expo- 
nentially growing contributions exp{±icj^t}. In the limit 
^(u>l )t ^> 1, the d k- integrals are dominated by the 
fastest-growing modes and can thus be approximated via 
the saddle-point method. Since the initial state (param- 
agnetic phase) is supposed to be globally isotropic such 
as the thermal ensemble or the ground state of |(2|), the 
dominant contribution to the two-point magnetization 
correlation function is calculated to be 



(F a (r)F b (r')) = C^6 ab 



exp{g cr t} J„(k max \r - r'\) 
Vi {k max \r - r'\) v 



(3) 



where J v is the Bessel function of the first kind with 
the index v = N/2 — 1. (The two-dimensional case was 
also considered in 0, Gi|.) Note that the above result is 
universal: Apart from u, q cl - and fc max , all information 
enters the pre-factor C™ ax only, which depends on q out , 
sweep rate, and temperature etc. 

Now let us study the creation of topological de- 
fects, i.e., 50(A r )-hedgehogs, by the symmetry-breaking 
quench. Their characteristic size, i.e., the extent of the 
paramagnetic core, can be determined by comparing the 
kinetic term oc V 2 with the C2-contribution in |l]) which 
yields £ s = \ j \J2\c2\Qm = \/^mq CT , i.e., the spin heal- 
ing length. In order to deal with well-separated and 
therefore well-defined topological defects, we assume that 
the dominant wave-length A max = 2ir/k max (determin- 
ing their typical distance, see below) is much larger than 
the spin healing length £ s , which amounts to requiring 
< <7cr — 2g out -C g C r- This condition is reasonably 
well satisfied in the experiment j3] where £, s f» 2.4 /im 
and A max w 16 /im. Within the saddle-point approxi- 
mation (i.e., focusing on the dominant modes), the term 
—V 2 — > fcm ax in the equations of motion resulting from 
([2]) can be neglected compared to g cr S> fc^ ax /m, and we 
may approximate F s» q^F/4. Consequently, at each 
spatial position, the magnetization F behaves as the co- 
ordinate of a harmonic oscillator (in N dimensions), be- 
coming inverted upon crossing the transition, cf. Fig. [TJ 
Splitting F up into its modulus F and unit vector of 
direction n via F = Fn, we find that F grows exponen- 
tially whereas the dynamics of n freezes: In view of the 
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conserved "angular momentum" L — F x F, we find that 
\h\ = \L\/F 2 (x F§/F 2 (t) decreases rapidly. Ergo, af- 
ter a short time (of order 1 / g cr ) , the magnetization F at 
each spatial position r grows exponentially in some given 
direction n(r), i.e., the system rolls down the parabolic 
potential hill whereby the direction of descent does not 
change anymore, cf. Fig. [TJ The correlations between the 
frozen directions n at different spatial positions are gov- 
erned by the initial state and determine the seeds for the 
creation of topological defects - i.e., objects which cannot 
be deformed to a constant n in a smooth way fl6| . 

The topological defects we are considering are SO (TV) 
(anti-)hedgehogs in TV spatial dimensions with the typi- 
cal order parameter distribution n = ±r/r and can be 
characterized by a non- vanishing winding number, i.e., a 
topological invariant or "charge" belonging to the homo- 
topy group ttn-i(Sn-i) = Z. The winding number is 
calculated by an integral over a TV — 1-dimensional hy- 
persurface enclosing a TV-dimensional volume [17]: 



F . F aQ~i... 



(4) 



Here e denote the Levi-Civita symbols and Sn-i = 
2tt n / 2 /T(N/2) the surface of the unit sphere in TV di- 
mensions. Since the winding number operator counts the 
difference between hedgehogs and anti-hedgehogs, its ex- 
pectation value vanishes, (9T) = 0, but its variance (Or) 
as a function of the enclosed volume yields the desired 
spectrum of net magnetization fluctuations. 

In order to calculate (Or), we make an additional ap- 
proximation: For a TV-dimensional harmonic oscillator 
F sa q 2 r F/A, the ground-state probability distribution 
p(F) cx F N ~ X exp{— F 2 } of the amplitude F is peaked 
around a (for TV > 1) non-zero value and the relative 
width of this peak decreases with increasing TV. There- 
fore, we approximate the operator F by an exponen- 
tially growing classical value F — > Fit) in all expectation 
values. This semiclassical approximation will become 
asymptotically exact in the limit TV \ oo, but we expect it 
to yield qualitatively correct results also in lower dimen- 
sions down to TV = 2, since p(F) is discernibly peaked 
even for TV = 2 and the topological defects are mainly 
determined by the (quantum) fluctuations of n (and not 
of F). For example, calculating [analogously to ([3])] the 
contribution of the fastest growing modes to the correla- 
tor (F 2 (t, r)F 2 (t, r')) = F A (t) + 2{F(t, r) ■ F(t, r')) 2 /N, 
we see that the accuracy of the semiclassical approxima- 
tion F — > F(t) increases for large distances \r— r'\ and/or 
large TV. Note that all the linearized F a possess indepen- 
dent initial ground states and commute with each other. 
Therefore, the n a do also commute with each other [so 
that operator ordering in Q is not an issue] and with 
F, which ultimately makes possible the semiclassical ap- 
proximation of F — > F . 

The above approach enables us to calculate (Or) in ar- 
bitrary dimensions TV > 2. After inserting |(4]) and using 
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FIG. 2: Scaling of the variance of the winding number (yt 2 ) 
with the sample size R. A logarithmic plot of (9r)/«; is used 
in order to visualize the asymptotic kItik scaling. 



the semiclassical approximation F(t,r) w F(t)n(r), the 
expectation value can be decomposed into a product of 
TV two-point correlators J3]) . Let us exemplify this proce- 
dure for two spatial dimensions: The transverse magne- 
tization can be decomposed via F± = F x + iF y = Fe 1 ® 
into its modulus F and phase $ which commute and are 
both self-adjoint ji^l- The winding number (@]) reads 21 1 
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Using Eq. (J3j) , replacing F — > F, and choosing the con- 
tour £ as a circle with radius R, we get [22] 



(0T 2 (i?)) 



J 2 (kx) . 



(6) 



with k — 2i?fc max , cf. Fig. [2] As anticipated, the charac- 
teristic length scale (e.g., typical distance between vor- 
tices) is set by the dominant wavelength A max (instead of 
the spin healing length £ s , for example) . 

The analytical expression ([6]) enables us to study the 
large-scale spectrum of the winding number, i.e., the dif- 
ference between the number of vortices and anti-vortices 
01 = 01+ — 0I_ . Two models for the scalingbehavior of 01 
are commonly discussed in the literature [2, EB| ■ Assum- 
ing a random distribution of vortices and anti-vortices 
(random vortex gas model), the typical difference |0T| 
scales with yJ%l±, Because the total number 0T+ + 0T_ 
increases proportional to the area R 2 , this model pre- 
dicts a scaling of (m 2 ) ~ R 2 (TV = 2). Alternatively, 
the random phase walk model assumes a random walk of 
the phase <& along the circumference. The accumulated 
winding number |0t| scales with y/R in this situation, 
implying a scaling of (OT 2 ) ~ R. Of course, both mod- 
els cannot be the final truth (consider the deformation 
of the circle to a slim ellipse, for example), and it is not 
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obvious how to generalize the random phase walk model 
to higher dimensions - but one can compare the different 
predictions with our approach. For large k, the expansion 
Ji(kx t oo) oc cos(kx — 3n/4)/^/~Kx yields (cf. Fig. ^ 

{m 2 (n > 1)) - KhlK, (7) 

which is close to the scaling predicted by the random 
phase walk model - but still not quite in agreement due 
to the additional factor of In k. This extra factor only 
occurs for TV = 2 and arises from the slow decay of the 
two-point correlator |(3]) at large distances \r — r'J -1 / 2 , 
which deviates from a local random phase walk [23J. Al- 
though the sample size in Q is probably not sufficient for 
extracting the factor In k, this prediction could be tested 
in future experiments. 

Let us study intermediate K values and compare our 
results with the experiment [?| by placing a circle with 
a radius of R = 5 /im into the center of the condensate 
cloud. Since R is significantly bigger than the spin heal- 
ing length £ s ss 2.4 /im but still sufficiently below the 
planar Thomas-Fermi radii of the cloud (12.8 /im and 
167 /im), the assumptions entering our derivation (such 
as homogeneity) should provide a reasonably good ap- 
proximation. For R = 5 /im corresponding to k ~ 4, 
we obtain (Or) ~ 1/3, i.e., the probability of finding an 
(anti-)vortex inside the circle is around 1/3. In view of 
the cigar shape of the cloud, we may place several circles 
in the central region of the condensate such that their 
mutual distance exceeds the correlation length (typical 
domain size) of A max /2 « 8 /im. Since these circles can 
be considered nearly independent, we would expect at 
least one vortex in the cloud - in reasonable agreement 
with 13. In contrast, the defect density is calculated in 
Ref. [19] to be fc^ lax /(47r), treating the magnetization as 
a Gaussian stochastic field. Therefore, a circle with a ra- 
dius of R = 2/fc max w 5/im should typically contain one 
defect. Although the discrepancy (1/3 vs. one) is not 
huge, it already suggests that our results are not quite 
compatible with Ref. and thus experiment should 
judge which approach yields the better description. Al- 
ternatively, numerical simulations (supplemented by suit- 
able assumptions about the initial noise spectrum) pro- 
vide a complementary approach [2J|. 

In conclusion, for the example of a spinor Bose gas, 
we studied the creation of topological defects during the 
quench from the paramagnetic to the S'0(A r )-symmetry 
breaking ferromagnetic phase. For an arbitrary number 
N of spatial dimensions, we derived a universal expres- 
sion for the winding number (which only depends on N 
and the enclosed volume in units of A^ ax ) and found a 
logarithmic correction to the scaling law of the random 
phase walk model for N = 2. 
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